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Abstract. We consider Gaussian states of fermionic systems and study the action 
of the partial transposition on the density matrix. It is shown that, with a suitable 
choice of basis, these states are transformed into a linear combination of two Gaussian 
operators that are uniquely defined in terms of the covariance matrix of the original 
state. In case of a reflection symmetric geometry, this result can be used to efficiently 
calculate a lower bound for a well-known entanglement measure, the logarithmic 
negativity. Furthermore, exact expressions can be derived for traces involving integer 
powers of the partial transpose. The method can also be applied to the quantum Ising 
chain and the results show perfect agreement with the predictions of conformal held 
theory. 
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Entanglement plays a key role in the study of quantum many-body systems BIZI. 
Considering a pure state of a composite system, a simple measure of the entanglement 
between two complementary parts is given by the von Neumann (or entanglement) 
entropy. Particularly interesting is the case of pure ground states where, in great 
generality, an area law for the entanglement emerges [3]. The most well established 
exceptions are one-dimensional quantum chains at criticality, where the entanglement 
entropy shows a universal logarithmic scaling [3] which can be fully understood with the 
help of conformal field theory (CFT) [2], The predictions of CFT have been confirmed 
on a variety of lattice models, among which a distinguished role is played by free-particle 
Hamiltonians. The ground states of these systems are given by bosonic/fermionic 
Gaussian states where the full entanglement spectrum is easily accessible |6J. 

The characterization of entanglement for mixed states is, however, far less obvious 
since, in contrast to the pure-state scenario, there is no unique way of defining a well- 
behaved measure. Among the numerous proposals for entanglement measures [7], a 
large family is based on a convex-roof extension of the von Neumann entropy. The 
drawback of these constructions is that they are essentially uncomputable already for 
systems of relatively small size. A viable alternative is based on an entirely different 
approach, making use of a special property of the partial transposition. Namely, the 
spectrum of the partial transpose of a density matrix may contain negative eigenvalues, 
only if the state is entangled In turn, a measure called logarithmic negativity [TO] 

can be introduced, which quantifies how much the partial transpose of a state fails to 
be positive, and can be shown to fulfill all the requirements of an entanglement measure 

HU- 

Although being a computable measure, the evaluation of logarithmic negativity 
might still pose a significant challenge in extended quantum systems. A notable 
exception is the case of bosonic systems, where the effect of partial transposition is 
equivalent to a partial time-reversal of the momenta in the corresponding subsystem H2t 
Furthermore, the partial transpose of bosonic Gaussian states remains to be Gaussian 
and, in turn, one has a simple formula to compute the logarithmic negativity via the 
covariance matrix [13] . Remarkably, the analogue statement does not hold for fermionic 
Gaussian states. 

The early studies of logarithmic negativity in lattice systems were conducted for 
the harmonic oscillator chain m ei m cm na m using the covariance matrix 
technique, and for spin chain models mm via density matrix renormalization group 
calculations. In contrast, exact analytical results were found only for a few simple 
spin models mm- A renewed interest in the problem was triggered recently, after a 
systematic approach within CFT was introduced [23]. This method could be applied to 
calculate the entanglement negativity for various geometries in ground [2T1I25] or thermal 
states m m of one-dimensional systems, as well as in out-of-equilibrium situations 
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Even though the predictions of CFT can be routinely tested on harmonic chains, 
calculating the logarithmic negativity in fermionic or spin systems remains to be 
more difficult. Recent studies employed a tensor-network representation of the partial 
transposition to calculate entanglement negativity for the transverse Ising chain [31] . 
Alternatively, Monte Carlo techniques were applied to calculate higher moments of the 
partial transpose [32, S3]. However, even for the simplest case of fermionic Gaussian 
states, a method which could compete with the computational simplicity of the bosonic 
case has so far been unknown. 

Here we show that, with a suitable choice of basis, the partial transpose of 
fermionic Gaussian states can be cast in a particularly simple form. Namely, it 
can be written as the linear combination of only two Gaussian operators, uniquely 
defined by corresponding covariance matrices which can be found explicitly. Under 
further assumption of a reflection symmetric geometry, this construction can be used 
to calculate a lower bound for the logarithmic negativity via the covariance matrix 
spectrum. For critical systems, the scaling behaviour of this bound shows remarkable 
similarities to that of the entanglement negativity. Furthermore, the higher moments of 
the partial transpose can be exactly evaluated through simple trace formulas, providing 
a way to test the universal CFT predictions on fermionic Gaussian states with minimal 
computational costs. 

In Section [2] we define fermionic Gaussian states and introduce the specific models 
in consideration. The partial transposition transformation for fermions is discussed in 
Sec. |3l focusing on a particular choice of basis which leads directly to our main result. 
Sec. E] is devoted to the construction of a lower bound for the logarithmic negativity, 
and its numerical investigation for the quantum Ising chain. Trace formulas for integer 
powers of the partial transpose of the reduced density matrix are presented in Sec. El 
The paper concludes in Sec. El with a short discussion of the results and their possible 
extensions. Various details of analytical calculations are included in three Appendices. 

2. Model and definitions 

We consider quantum systems associated to free-fermion Hamiltonians 


AT 



( 1 ) 


where the matrices A and B are Hermitian and antisymmetric, respectively. 
The fermionic creation/annihilation operators, and c m , satisfy the canonical 
anticommutation relations = S mn . For our purposes, it will sometimes be 

more convenient to work with Majorana fermions defined as 



( 2 ) 


satisfying the relations {a*,, a{\ = 2 5ki- In terms of Majorana operators, the Hamiltonian 


of Eq. ([T]) with real A and B can be rewritten as H = i^f 

^2m,2n— 1 -^ 2 n—1,2m 4 (dmn"fHmn) aild T2 m ,2n -^ 2 m—l, 2 n —1 0- Tile product of all 
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Majorana operators define the parity operator, P = i N n«=i a «) which plays an 
important role in fermionic systems. According to the parity superselection rule, only 
density matrices that commute with P correspond to physical states [Ml EM. 

The states we are going to study in this paper are the so-called Gaussian states. 
These describe the ground and Gibbs states of quadratic Hamiltonians, and play a 
prominent role in quantum information theory [37, 38] [39] . A state p is Gaussian if it 
can be written as 


p = — exp(^ W k ia k ai/A ), 

k,l 


(3) 


where IT is a purely imaginary antisymmetric matrix (with the possibility of \ W k i\ —> oo 
allowed) and Z is the normalization factor. A Gaussian state can also be characterized 
uniquely by its covariance matrix T w = ([a^, )/2 via 

W 

tanh — = T. (4) 

Using the covariance matrix, one can express the expectation value of any Majorana 
monomial through the Wick expansion 


i 

Tr(p 

^ j n\^ j n2 • • • E sgn (7r) n ^ n Tt{2k-T), n n(2K) > 

7T fc= 1 


(5) 


where all the indices are different, the sum runs over all pairings n, and sgn (n) denotes 
the sign of it [t|. Let us note, that similarly to Gaussian states, one can introduce general 
Gaussian operators which are also defined through Eq. (|3]) , however, without requiring 
that the spectrum of W is real. The Wick expansion, i.e. Eq. (ED, holds for these 
operators, as well. 

We will also study spin chain models that are related to free-fermion Hamiltonians 
through the Jordan-Wigner transformation [4U] 

m— 1 m— 1 

a 2 m-l = n a k a m’ °2m = J(6) 

k =1 k =1 


where (with a — x,y, z ) denote the Pauli matrices acting on site m. The prototypical 
example is the transverse field Ising (Tl) chain, 


Htt — 


E( 


< An ( An+1 


+ , 


(7) 


where a chain of length N with open boundary conditions is considered. Applying the 
Jordan-Wigner transformation, the Tl Hamiltonian ([7]) takes the form of Eq. (JTD with 


-A-mn cy (&m,n— 1 T $m,n+ 1) hS n 


Bmn 1 $m,n-T l)- 


( 8 ) 


f A pairing n over a set (1,2,..., 21} is a permutation of the 21 elements which satisfies 7r(2fc — 1) < 
ir(2k) and n(2k — 1) < n(2k + 1) for any 1 < k < L Any pairing can be decomposed into an M 
number of transpositions (simple exchange of only two indices), and the sign of the pairing is defined 
as sgn (7r) = (—1) M . 
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Such a quadratic Hamiltonian can be diagonalized by a canonical transformation, 
N 1 

Vk = ^2 9 [^( m ) ° 2 m-i - i'f’kim) a 2m \ , 

ra=l 

and brought into the standard form 

N 

H = ^2 ^kvivk + const. 

fc=i 


( 9 ) 


( 10 ) 


The spectrum A k and the vectors (fikif’k in Eq. (JH]) follow from the eigenvalue equations 

(A-B)(A + B)(f> k =Al(i> k , (11) 

(A + B)(A-B)if k = Aj^ fc . (12) 


With the full solution of the problem at hand, one can directly write down the covariance 
matrix for a Gibbs state, with inverse temperature /3, of the open TI chain as 


fin nia • ■ • niAr 

1I21 n 22 


n^i . n NN 



(13) 


with matrix elements given by 


Qmn 



ifk{rn)(j)k(n ) tanh 


2 


(14) 


k 

In our studies we will also be interested in the periodic TI chain, given by a 
Hamiltonian as in Eq. (J7]) with the sum running up to L and boundary condition 
a x L+l = af. Note that, for clear distinction, we will use L instead of N for the length of 
the periodic chain. It is well known, that its ground state is the same as for the fermionic 
model with matrices as in (jHJ) but with antiperiodic boundary conditions. The solutions 
of the system (HID and (fT2|) are plane waves, 0fc(m) ~ e Wkrn and r f k ( rn ) ~ e* ei: e* PA:m , with 
the Bogoliubov angles and eigenvalues given by 


h — _ 

e l6k = —r-, Afc = sjl + h 2 - 2/icos pk, (15) 

and the allowed values of the momenta are p k = (2k— l)7r/ L with k = —L/ 2 + 1,..., L/2. 
One can also work directly in the thermodynamic limit, L —» 00 , where the momenta 
become continuous and the sum in Eq. (1141) for the matrix elements g mn is replaced with 
an integral. 


3. Partial transpose for free fermions 

As discussed in the Introduction, the partial transposition plays an important role in 
quantum information theory. In the context of entanglement theory, it was first studied 
for qubit and qudit systems M, but later also bosonic H21 HD 021 @3] and fermionic 
models [3S EH El 05] were investigated. An important result coming from these 
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studies was that the partial transpose of a bosonic Gaussian state is again a Gaussian 
operator; this simplifies the calculation of the negativity [TBJ HSJ ■ The analogous result 
for fermionic Gaussian states does not hold, which can already be demonstrated by 
2-site systems, see Appendix A 


This section is devoted to the derivation of a weaker, but still useful, result for the 
fermionic case. After recalling the notion of the partial transpose for spin systems and 
the corresponding definition for fermions in Section IHTT1 we show in Section 1X21 that the 
partial transpose of a Gaussian state (in a particular basis) can always be decomposed 
as a sum of two Gaussian operators. This decomposition lies at the heart of all the 
results in the further sections. 


3.1. Definition of the partial transpose 


B Aj B A 2 B B A! A 2 B 


(a) (b) 

Figure 1. Two possible partitioning of a spin/fermionic chain into subsystems A\ 1 
A 2 , and B, as described in the text. 


Consider a general tripartition of a chain of qubits into disjoint sets A 1 , A 2 and 
B , e.g. as in Fig. [0 Let p denote the density matrix of the whole composite system. 
Defining A = Ai U A 2 , the reduced density matrix (RDM) of subsystem A is given by 
Pa = Tr B p. The partial transpose of the RDM, p T f , with respect to the subsystem A 2 
is defined by its matrix elements as 


( e i e i \P/ 


\ e k e ' 


>) = {e^e^p^e ?), (16) 

where {|e J - 1 ' ) )} and {| e ( p )} denote complete bases on the Hilbert spaces "Hi and l-i 2 
pertaining to the subsets Ai and A 2 . The definition of pj is basis dependent. However, 
one can easily characterize the set of transpositions on the operators acting on 'H 2 as 
those non-degenerate linear transformations that satisfy 


n{M x M 2 ) = 1Z(M 2 )1Z(Mi) , 


(17) 


for any two operators M\ and M 2 acting on l-i 2 . Since any two partial transpositions 
can be connected by a unitary conjugation, the eigenvalues of p T Jf are independent of 
the choice of basis. Moreover, it was shown that the partial transpose of the density 
matrix can only have negative eigenvalues if the corresponding state is entangled PEj- 
In a similar way, one can define the partial transpose for fermionic states. Consider a 
tripartition of a system with N fermionic modes, e.g. as in Fig. |T]. Let {mi, m 2 ,..., m 2 k} 
and {ni,n 2 ,... ,n 2 ij denote the indices of the Majorana operators belonging to the 
subsystems Ai and A 2 , respectively. Let us introduce the notation a x = 1 and a l x = a x . 
A general fermionic state on A = Ai U A 2 can be written as 

p A = Y^ w *,z a ^ • • • a™ , 

k,t 


( 18 ) 
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where the variables k = (k\, ... K 2 k) and r— (r\,, T 21 ) in the summation run over all 
bit-strings of length 2k and 2£, respectively. Note that since physical fermionic states 
must commute with the parity operator, as discussed in Section [2j one has Wk, z =0 when 

E 2 k , • in 

i=\ Ki + E ;=1 T j 1S odd - 

The partial transpose of pa is simply a transformation that leaves the A\ Majorana 
modes invariant and acts as a transposition on the operators built up from modes of 
A 2 , i.e. 


Pa = V ... <C” X( < ... <“). (19) 

where 1Z satisfies Eq. CGJ- Since also in the fermionic case all the transpositions are 
connected by a unitary conjugation, the eigenvalues of pj will be independent of which 
1Z we choose. It will be useful to consider the following particular transposition which 
is defined by 


1Z{ 


u ni * * * u n 2 e) 


(-1 y iL) a T n \ ... a™, where /(r) 


0 if Irl mod 4 G {0, 3}, 

i-i Vi 2Q 

1 if |r| mod 4 e {1, 2}, 


where |r| = Ei=i r *- dl other words, a Majorana monomial is mapped by 1Z to itself if 
it is of length 4n or 4n + 3, and otherwise it is multiplied by a —1 sign. Note that a very 
similar transposition in fermionic systems was already considered in Ref. 03. Although 
the definition of entanglement in fermionic systems is somewhat different from the case 
of spin systems, it has been proven that p T jf can only have negative eigenvalues if pa is 
entangled also in the fermionic case }34l 135], 

Finally, let us shortly discuss the connection between reduced density matrices of 
fermionic and spin models that are connected by the Jordan-Wigner transformation. As 
this transformation is non-local, it has been shown that the reduced density matrices 
corresponding to a region A\ U A 2 in a spin chain model and its fermionic counterpart 
are usually not equivalent (not isospectral), unless A\ and A 2 are adjacent intervals, as 
depicted in Fig. [T](b) [4H1 EE]- Since the same holds also for the transposed density ma¬ 
trices, when treating spin models we will only consider the adjacent interval geometry. 
For the case of fermionic systems, our results are valid for arbitrary geometries. 


3.2. The Gaussian case 


Consider a Gaussian state pa on the system A 



A 1 U A 2 with a covariance matrix 


( 21 ) 


where T 11 and T 22 denote the reduced covariance matrices of subsystems A\ and A 2 , 
respectively; while T 12 and T 21 contain the expectation values of mixed quadratic 
terms. Let P 2 be the parity operator on subsystem A 2 , and define the operators 
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p + = KPA+P 2 PAP 2 ) and p_ = \{pA—P 2 pAp 2 )- By definition, we have that pa = P++P-- 
Using the notation of Section IHTT1 p + and p_ can be expanded as 

P+= w ^ a ™i • • • a %2k a n\ ■ ■ ■ a nL » 

|r| even 

p-= J2 w ^- < • • • • • • ° r z > 

K , T _ 

|r| odd 



where the coefficients w^,t can be obtained from R 4 using the Wick rule, Eq. (J5J). By 
linearity of the partial transpose, pj = p + 2 + p Tz follows, and p± can be obtained using 

Eq. ([20]): 


„t 2 _ 


p+ = 5 ^ (-1 • • • « T2£ 


«,r 

|r| even 


p- 2 = ( _1 ) 


(IXI —1)/2 K! 

w k,t “mj 


1 


q t 1 

U m 2k a n 1 


n 2 ^ 1 


7 t 2 i 
l n 2t 


(23) 


|r| odd 

Let us introduce the generalized Gaussian operators 0 + and CL, with covariance 
matrices 


rv = 


r 11 fr 12 
ir 21 -r 22 


r_ = 


r 11 -*r 12 

-fr 21 -r 22 


and consider the Majorana monomial expansion of these operators, 

0± = °t: 


n 1 n ZK n x n 
l,r a mi * * • a m 2k a ni * * * a n 2 e 


(24) 


(25) 


«,r 


Since 0 + and CL are Gaussian operators, one can again obtain of T from r± using 
Eq. ([5]). Connecting the Wick-expansion form of w^ z with that of of T , using the relation 
between IL and r±, one can deduce that 


Otr = 


f ±i(-l) (l - 1 1)/2 MU,r 

when r 

\ {-IpPw^r 

when r 


(26) 


Comparing this with Eq. (125]) it immediately follows that p + 2 = \{0 + + 0_) and 
pj 2 = f(0- — 0 + ). Thus, we obtain the decomposition 


P T 1 


1 — i 


CL 


1 + i 


CL, 


2 2 

which is, from a conceptual point of view, the main result of the paper. 


(27) 


4. Partial transpose and logarithmic negativity 

In the previous section, we have shown that the partial transpose of a Gaussian RDM 
can be written as a linear combination of only two Gaussian operators, which is the 
simplest possible form for a non-Gaussian operator. However, since 0 + and CL do not 
commute in general, Eq. (127|) can not be rewritten for the eigenvalues, and thus one 
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does not have a direct access to the full spectrum of pj. Nevertheless, there are a 
number of important properties which can be deduced by a direct investigation of the 
covariance matrices T±. A particularly interesting quantity that we will study is the 
logarithmic negativity uni, which can be used as a measure of entanglement. 

In the following we thoroughly investigate three special cases. First, we consider the 
partial transpose for bipartite pure states which, although the results being well-known, 
turns out to be very instructive in understanding the implications of the decomposition 
in Eq. (127[) . We proceed with the study of thermal mixed states in a reflection symmetric 
bipartite geometry, which allows us to define and calculate a lower bound for the 
logarithmic negativity. Finally, we report our findings for a genuine tripartite geometry. 
In each of the following subsections, the validity of formula (127)) was checked against 
exact numerical calculations for the TI chain, Eq. (17)1 . with a small number of spins. 


4-1. Bipartite pure states 


We first consider the simplest case with B = 0 and a pure state on A = A 1 U A 2 given 
by p = | 0 ) (0|. Since for any pure Gaussian fermionic state T 2 = 11 is satisfied, it follows 
that [r+,r_]=0, which implies that 0 + and 0_ commute as well. Furthermore, since 
their eigenvalues are connected by a complex conjugation, the spectrum of p Ta is simply 
given as the sum of the real and imaginary parts of the 0 + eigenvalues. Now, since 0 + 
is also Gaussian, its spectrum is uniquely determined through the eigenvalues of T + . To 
obtain them, we first assume without loss of generality that |Ai| < \A 2 \. Furthermore, 
for notational simplicity we also assume that |A| is even, however, the results for the odd 
case follow trivially. Let ±/ifc denote the eigenvalues of the reduced covariance matrix 
T 11 with k — 1 ,..., | Ax| and p, k > 0. As shown in Appendix B the eigenvalues ±z/j^ of 
r + can then be given as 


v k = 


p/c ± 1 — Pfc when k = 1 ,..., |Ai|, 

1 when k = |Ai|+l,..., |A|/2 . 

The canonical diagonalised form of the Gaussian operator 0 + reads 

0 = tt tt l + 

k=1 cr/ c =d= 


(28) 


(29) 


where bf are Majorana operators obtained from aj via the orthogonal transformation 
which diagonalizes T + . The eigenvalues of 0 + can be obtained according to the following 
rules. First, we shall consider the various combinations of the conjugate eigenvalue pairs 
for each k = 1 ,..., |Ai| as 


OJu = 


1 + <y k v£ 1 + <r' k v k 


(30) 


with = ± and a' k = ±. Using Eq. 


/ .++ _ 
— 


1 + pfc 


OJi. = 


this yields 




= -U, = 


i-/4, 


2 


2 


(31) 
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where we used the property vffu k =1. The nonzero eigenvalues of 0 + can then be 
written down as 

sw=rk K = n 1 - ± y b - n 

fc =l a k=v' k Vk = -&' k 

where g_ and of are the signature arrays corresponding to the eigenvalue. Note that the 
additional off — 1 in Eq. (125j) lead to further eigenvalues of 0 + that are all equal to 
zero. 

The products in Eq. (152]) are either real or purely imaginary and the eigenvalues of 
p T2 thus follow as Re£\CT/ or respectively. It is instructive, however, to derive 

the same spectrum using the Schmidt decomposition of \<f). Dividing the Hilbert space 
Ti — Tii <S> %2 into two parts [jj, one has 

W=£Al#>l*f>, (33) 

i 

with the RDM eigenvalues A* of subsystem A\. Clearly, the state is supported on a 
smaller Hilbert space "Hi <E> "Hi and is invariant under the action of a flip operation 
defined by |0 |)|$■) —> \f>))\ 0f). The partial transpose of p, 

p T2 = (i0K0i) T2 = ^\m)\ ® ( 34 ) 

i,3 

commutes also with the flip operator. Furthermore, it is easy to check that the 
eigenvalues and vectors are 

P T2 I0])l0i) = M 0])l0i), P T2 \^j) = ±\A^|0$), (35) 

where we introduced the notation 

10$) = ^=G0i)l0?) ± 10]) 10?)) * + j- (36) 

Note that all the positive (negative) eigenvalues correspond to even (odd) eigenvectors 
with respect to the flip operation. Moreover, since p is Gaussian, one can immediately 
write down the products of eigenvalues as 

vCV = n l±£Wf. U £EK, (37) 

a k =cr ' k a k = ~ cr ' k 

where the signature a has again components cr fc = ±. Comparing Eq. (157)) to Eq. fl32|) . 
one indeed recognizes up to the factors of i. 

Owing to the simple product structure of the eigenvalues in Eq. (1571) and, in 
particular, to the fact that all the negative eigenvalues are located in the odd subspace, 
one has 

Tr \p T f = 1 - 2Tr 0 / 2 = (l + yjl - , (38) 

§ When considering a tensor product structure, we always refer to a partition of the spin chain system 
constructed through the Jordan-Wigner transformation. 
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where we introduced the notation Tr e / 0 for traces taken over the even/odd subspace. 
Thus, the pure-state logarithmic negativity is given by the simple formula 

e = InTr l/’l = £ In (l + fl - ff) ■ (39) 

k ' ' 

Considering also the expression of the Renyi entropy for fermionic Gaussian states, 

1 ~ hk \ 

2 ) 

the well-known equality £ = S' 1/2 for pure states can be confirmed directly. 


(40) 


S„ = 


n 


S ln 


1 + hk 


+ 


f.2. Thermal states in a bipartite geometry 


The simplicity of the pure-state scenario relies essentially on the property of the Schmidt 
decomposition, which is automatically symmetric under the flip operation defined 
previously. Due to this, the partial transposed state is block diagonal wrt the splitting 
into even and odd subspaces. Such a structure is missing for general bipartite mixed 
states, unless the system has a flip-type symmetry a priori. In this respect, a natural 
scenario would be to consider intervals of equal length |Al| = \A 2 \ = N/2 and states 
that are reflection symmetric. To analyse such a situation, we shall consider Gibbs 
states of the open TI chain, Eq. (j7|l . with a symmetric bipartitioning, these being the 
simplest mixed Gaussian states where we hope to get further insight into the structure 
of the partial transpose. 

Considering a covariance matrix of the form the spectrum of T + is invariant 
with respect to a sign change and complex conjugation, hence the eigenvalues can be 
collected into two families of quadruplets 


{z k , z* k , -z k , -z* k } , k G (I) {iu k , -iv k , -iu k , iv k } , k G (II) (41) 


where in family (I) we choose Rez*, > 0 and Im z k > 0 , whereas u k > v k in family 
(II). Note that the eigenvalues in the second family are purely imaginary and thus their 
complex conjugate are automatically contained in the spectrum of a skew-symmetric 
matrix, hence u k 7 ^ v k . Although one could, in general, have an arbitrary number of 
type (II) quadruplets, from the numerics we observe that in the Ising case they are either 
absent or a single one appears. Moreover, this only happens in the symmetry-broken 
phase, i.e., when h < 1 in Eq. Q. 

Analogously to the pure case in Eq. (l30lh we first assign the factors 



k 


1 + \zk I 2 + cr k 2Rez k 
4 

1 - \z k \ 2 + a k 2ilmz k 
4 

1 + u k v k + a k i(u - v ) 
4 

1 - u k v k + a k i{u + v) 
4 


= cr'k 


°k = - 0 *. 


fce(l) 


Vk = <j' k 
Vk = -Ofc 


k G (II) 


( 42 ) 
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within each quadruplet, and the eigenvalues £lag> of 0 + are again given in the factorized 
form of Eq. (1521) . Although the spectrum of the operator 0_ is identical to that of 0 + , 
they do not commute in general and thus one has no direct access to the eigenvalues 
of p Tz . Nevertheless, the information about the even/odd parity of the eigenvectors is 
retained. In fact, the reflection operator R , which defines the even/odd subspaces in 
our case, acts on the spin operators as Ra^Rf = cr^_ n (with a = x,y, z), implying the 
action Rc\R t = Pcjy_ n on the creation operators, where P is the parity operator. Using 
this, it follows that the sign factor associated to an eigenvector of 0 + reads 


Sg_a' 11 ll S<T fc <T U Sa k a 'k 



if a k = a' k , 
if a k = -a' k , 


(43) 


k k 

which can also be verified by considering the pure state limit. The parity of the 0_ 
eigenvectors are simply obtained through the factors s* 




With the knowledge of S arT ', we are now able to carry out signed traces of the form 


Tr c 0 + - Tr 0 0+ = SaA*’ ■ (44) 

a,a f 

Note that the terms in the sum of Eq. (145)) completely factorize in the quadruplet 
index k. Furthermore, using the fact that Tr c / 0 0_ = (Tr 0 / 0 0 + )* and thus Tr e / 0 p T2 = 
Re Tr e / 0 0 + + Im Tr e / 0 0 +1 a simple calculation leads to 

1 - 2Tr 0 p T * = Tr e p Ti - Tr= [[ l + kl 2 + 2Imz fc j-j 1 + u t v k + u k + (4g) 

fee (I) fce(ii) 

Finally, we define the quantity 


£ 0 = lnmax(l — 2Tr 0 p T2 ,1), 


(46) 


which clearly gives a lower bound for the logarithmic negativity, £ > £ 0 , with strict 
equality if all the negative eigenvalues reside in the odd sector and their number is equal 
to the dimension of that subspace. This is true for pure states and one expects it to be 
valid for thermal states in a finite regime of temperatures above the ground state. For 
general bipartite states, however, the dimension of the negative subspace can be much 
larger [50j [5Tj. 

To test the bound £ 0 against the exact value of the logarithmic negativity, we 
considered small TI chains with iV < 10 and obtained £ via exact diagonalisation of 
p T2 . This is shown on Fig. [21 as a function of the temperature (left) as well as of the 
magnetic field (right). We find that, for low enough temperatures, £ 0 indeed exactly 
coincides with £. For larger temperatures, however, some of the negative eigenvalues in 
the odd sector become positive or, vice versa, even eigenvectors could attain negative 
eigenvalues, with both of these processes increasing the difference £ — £ a . Nevertheless, 
for the small system sizes considered, it appears that £ 0 gives a very good approximation 
up to temperatures T fa 1/N, above which it starts to deviate significantly. 
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Figure 2. Logarithmic negativity £ (symbols) vs. £ 0 (lines) between two halves of 
a small Ising chain with Hamiltonian 0 in a Gibbs state. Left: as a function of the 
inverse temperature /3 with h = 1 and various number of spins N. The curves £ and 
£ 0 start to deviate around j3 rts N. Right: as a function of h for various /3 and N = 8. 
Minor deviations between £ and £ a are only visible for /3 = 10. 


f.3. Ground states in a tripartite geometry 


So far we have only considered bipartite geometries. Another interesting setting we are 
able to deal with is a pure state with the symmetric tripartite geometry, depicted on 
Fig. QJb). In this case, the reduced state pa after tracing out the sites of B is a mixed 
Gaussian state, associated to the reduced covariance matrix with indices running 
over sites in A [52j. 

The logarithmic negativity for the tripartite case can be obtained with CFT 
methods pSj. For two intervals of the same size i embedded in a system of length 
L with periodic boundary conditions, the calculation yields m 


£((, L) = - In 


L sin 


+ const., 


(47) 


vr sin (^) 

with the central charge c and a non-universal constant. However, subtracting the value 


at i = L/A one obtains a universal scaling function 

e(z) = £(£, L ) — S[Lj 4, L) — - In [tan(7 r*)], 


(48) 


where z — i/L and we have set c = 1/2 corresponding to the TI chain. The formula 
was tested using tensor network methods for the calculation of the partial transpose, 
and a very good agreement was found [HI] . 

It is interesting to check the behaviour of the lower bound, defined in Eq. (j45|) . for 
the geometry at hand. In fact, since reflection symmetry is fulfilled, all the arguments 
of the previous section, leading to Eq. ([45]) , apply and £ 0 is given by the same formula 
in terms of the eigenvalues of r + . However, there is an important difference compared 
to the bipartite thermal case, which is apparent from the numerical investigation of 
small systems. Namely, the number of negative eigenvalues of p T f is always less then the 
dimension of the odd subspace and, moreover, some of the corresponding eigenvectors 
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are even. Thus, by tracing out the sites of B , the partial transpose pj cannot be 
smoothly deformed from the pure-state case and, consequently, one does not have a 
finite regime of parameters where the bound given by £ 0 is tight. 

In spite of the above findings, £ 0 shows a very interesting behaviour, which is 
demonstrated on Fig. [3j First of all, it shows a clear signature of the phase transition 
at h — 1, which can be seen when plotting £ 0 (L/4, L) against h on the left of Fig. [3j 
Furthermore, defining the quantity e Q analogously to Eq. (1451) . one finds an excellent 
data collapse when plotted against the variable z = £/L, see right of Fig. [3j The scaling 
function is found to be given by 


e 0 (z)=£ 0 (£,L)-£ 0 (L/4,L) = —]n 

16 


tan(7 tz) 


(49) 


2 COS 2 (7T2) 

Although the functional form of e Q (z ) was found by trial, one has an excellent match 
with the data without any fitting parameters involved. Interestingly, the prefactor of the 
logarithm is exactly the half of e(z) in Eq. (1481) . however, the argument is modified as 
well. We also performed a calculation directly in the thermodynamic limit L — > oo, and 
found £ 0 (£) = 1/16ln£ + const., which is perfectly consistent with the above findings. 
Furthermore, one could also consider the simple fermionic hopping chain (or XX chain in 
spin language), defined by B mn = 0 and A mn = y(<5 mjn _i + <5 m ,n-i-i)- In complete analogy 
with the result for the bipartite entanglement [53], one finds £ xx {2£) = 2£ff I [£f) for 
h — 1, and thus a doubled prefactor 1/8 with respect to the critical TI chain. Therefore, 
even though £ 0 does not approximate £ well, it shows exactly the same universal 
behaviour, suggesting it as an entanglement indicator which is extremely simple to 
calculate. 



Figure 3. The scaling of £ 0 (i. L) for two adjacent intervals of equal length l in the 
ground state of a periodic Ising chain with L sites and magnetic field h. Left: £ 0 {L/ 4, L) 
as a function of the magnetic field h for various L. The logarithmic divergence around 
h = 1 is clearly seen. Right: £ 0 (£,L) — £ 0 (L/4,i) at the critical point h = 1 against 
the scaling variable z. The solid line shows the scaling function in Eq. 
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Our result for the partial transpose, Eq. (I27|h can be further tested by looking at traces 
of integer powers of the partial transpose p T A which can also be carried out within CFT. 
Since the identity Tr(pJ ) 2 = Tr p 2 A holds for any density matrix, the simplest nontrivial 
quantity to check is the trace of the third power. For the geometry of the previous 
section, one finds the CFT result 


R 3 (£,L) = \nTr(p T J) 3 = --ln 


L 3 


Sill 


7T J 


7 ll 

T 


Sill 


27 tI 


+ const. (50) 


Similarly to Eq. (148 H . a universal scaling function can be defined as [M] 

r 3 (z) = R 3 {£, L) - R 3 (L/4, L) = ~ In [2 sin 2 ( 7 rz) sin(27r^)] , (51) 

which was already tested numerically for the critical TI chain [31, [32]. On the other 
hand, one could also consider two adjacent intervals of equal length £, embedded in 
an infinite chain which is thermalized at inverse temperature j3. Applying a simple 
conformal transformation, the corresponding CFT formula follows as 


R 3 (t,f3) = 


9 


In 


■/3 3 


7T 


sinli 


7t£ 

J 


sinh 


2t t£ 

~T 


const. 


(52) 


We now show how the above traces can be calculated with the covariance matrix 
formalism. Expanding the third power of p T f in Eq. (I27|) and taking the trace one 
arrives at 


Tr(/>5) S = + |Tr(0^0_), 


(53) 


where we have used that both of the traces on the right hand side are real. In order to 
evaluate them, one has to invoke the determinant formulas for the trace of products of 
Gaussian operators, which have already been considered in different contexts mm- 


The main steps of this calculation are summarized in Appendix C In turn, one finds 


Tr(ri s r = =Fj 



(54) 


where the sign of the first term depends on the spectrum of r + , see Eq. W- Namely, 
the + sign applies only if the quadruplet with purely imaginary eigenvalues appears 


see 


Appendix C for a more detailed discussion). Note that similar sign ambiguities 
One should also remark, that traces of higher powers can be 


also appeared in 

handled in a very similar way, however, the formulas become rather lengthy. 

The trace formula (l54|i can now be compared to the CFT predictions in (T5TT) and (1521) 
by inserting the corresponding covariance matrices E ± and evaluating the determinants. 
This is shown in Fig.[4]for the ground (left) and thermal states (right), respectively. The 
perfect agreement of the curves provides a highly nontrivial check of the CFT results. 
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Figure 4. Left: R 3 R, L) — R^{L /4, L) as function of z = Z/L for two adjacent intervals 
of length Z in the ground state of the critical TI chain of size L. The solid line shows 
the CFT formula (l5lT) . Right: R^{Z, (3) for adjacent intervals of size Z in a thermal 
state of the infinite chain (h = 1), with inverse temperature (3. The inset shows the 
rescaled data compared to the CFT prediction (15^1) on a horizontal log-scale. 


6. Discussion 

In conclusion, we have shown that the partial transpose of fermionic Gaussian states can 
be written as a linear combination of only two Gaussian operators, uniquely determined 
by associated covariance matrices T±. I 11 the presence of reflection symmetry, this 

particular form of the partial transpose allows us to carry out traces over the even/odd 
subspaces which, in turn, can be used to construct a lower bound to the logarithmic 
negativity. Furthermore, the trace of any integer power of p T ff can, in principle, be 
calculated as a sum of determinants, each of linear size 2\A\. 

There are several open questions left for future research. We did not consider 
in detail entanglement detection questions, e.g., providing temperature bounds for 
separability of fermion or spin systems in thermal equilibrium. It would be instructive 
to compare such results obtained from the negativity lower bound £ 0 with earlier studies 

[55 , "561 157] . 

Another natural extension of our work would be to consider non-adjacent intervals. 
For spin chains, however, it was shown that the spin RDM itself is already a linear 
combination of four fermionic RDMs [48] , Although our construction for the partial 
transpose could be carried over, it would further double the number of terms in the 
linear combination. Thus, the calculation of the traces for such a case is still realisable, 
but presumably more tedious. 

It would also be interesting to see whether the lower bound £ 0 could be attainable 
within the framework of CFT. This could lead to an analytical understanding of the 
scaling function for the critical tripartite case in Eq. (1491) and could shed light to the 
origin of the prefactor. In fact, one is tempted to guess that this is equal to one-half of 
the corresponding prefactor of the logarithmic negativity in a general CFT. From the 
free-fermion point of view, our analysis clearly suggests that certain asymptotic relations 











On the partial transpose of fermionic Gaussian states 


17 


between £ 0 and £ could hold in general. Finding a rigorous form of this relation would 
allow for a numerically feasible estimation of the entanglement negativity for fermionic 
systems. 

Finally, we point out that some specific classes of mixed Gaussian states exist 
which allow for an exact calculation of the logarithmic negativity using the methods 
introduced here. These are the states for which the relation [r + ,r_] = 0 is satisfied, 
an example being the isotropic Gaussian states [33], for which the covariance matrix 
satisfies T 2 = A 2 11 with some 0 < A < 1. The situation is then analogous to the pure- 
state case and the corresponding calculation can be generalized, which we leave for 
future studies. 
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Appendix A. The partial transpose of a 2-site RDM 


It is instructive to check how the method introduced in Section [3] works for the simplest 
case of two consecutive sites. The canonical form of the Gaussian RDM is given by 

n il + ZZ/fc6 2 fc-1^2fc / a , X 

- 2 -’ (A- 1 ) 

fc=l,2 


where zy, ; £ [0,1] and the Majorana operators bj are related to aj through an orthogonal 
transformation. For simplicity, we will consider only covariance matrices of the form 
of Eq. 013 p , and assume also reflection symmetry. These states are parametrized by a 
single angle 9 beside the covariance matrix eigenvalues 

Using the Jordan-Wigner representation of aj and working in the usual spin basis, 
the most general form of the partial transpose is 


P T2 = l + 


11 v\v 2 


—Z —Z | 

cr 1 a 2 + 


zq + v 2 


( cos 26 

V 0 


o 

sin 29 


sin 26 
0 


0 ^ 

+ " 1 ;" 2 

( 0 

0 0 

i \ 


4 

0 0 


cos 26 J 


W 

°) 


(A.2) 


Note that, besides the diagonals, all matrix elements vanish and are thus not shown. It 
is straightforward to obtain the four eigenvalues 


Al,2 — 

<^ 3,4 


! + v x v 2 ± [y 1 + zq) 2 COS 2 29 + (z /1 - z/ 2 ) 2 

[1 - V\V 2 ± (zq + v 2 ) sin 29} /4. 


/4, 


(A.3) 
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Using the parity operator P 2 = erf, one can also construct p ^ 2 = (p 72 ± -P 2 P 22 P 2 )/2 with 
matrix elements 


t 2 i , ^2 
=4 + — 


T, + ^2 




/ cos 26 



0 

\ 



zq + zq 


0 

0 




°"l°2 + 

4 


0 

0 



5 



V 0 



— cos 26 

) 


/ 0 

0 ^ 




( ° 


1\ 

0 

sin 26 

Ul 

— 

zq 

0 

0 


sin 26 

0 


4 


0 

0 


v° 

°) 




l 1 




(A.4) 


(A.5) 


The eigenvalues of the operator p+ + ip 72 then read 

^ 1 , 2 = 1 + zqz / 2 ± \/{yi + zq ) 2 cos 2 29 - [yi - z / 2 ) 2 /4, 
U2 3 4 = [1 — zqz / 2 ± i(yi + z/ 2 ) sin 29] /4. 


(A.6) 


Note that we have A 3i4 = Re Q 3)4 + Im Q 3i4 which simply follows from the fact that 
pff and p Ta commute on the corresponding subspace, including the odd eigenvector. 
Unfortunately, this property does not generalize to symmetrically bipartitioned intervals 
with more than two spins. 

We will now show that the Gaussian operator 0 + with covariance matrix T + has 
indeed eigenvalues given by Eq. (1A.6[) . The covariance matrix T for the Gaussian state 
(lA.ip and the associated T + have the form 


/ 0 

c 

0 s_ \ 


( 

0 

c 

0 

\ 

—c 

0 

s + 0 

, T + = i 


—c 

0 

is + 

0 

0 

-8+ 

0 c 



0 

—is .|_ 

0 

—c 

V - s - 

0 

—c 0 j 


V 

— 

0 

c 

0 ) 


T = i 


with the shorthand notation 
+ v 2 


c = 


cos 29, 


s± = 


+ V2 


sin 29 ± 


v\ - v 2 


(A.7) 


(A.8) 


The four eigenvalues iz/^ of T + can be computed with 


v 


± 




,S+ + S- 


(A.9) 


If the operator 0 + is Gaussian, its eigenvalues must have the form 

1 + „ + 1 + „ l-v+l-v 

2 2 2 2 

sy = 1 + v + 1 ~ v - i -. + i + ^ 

2 2 2 2 


(A.10) 


Substituting (IA.8D and (IA.9[) into flA.101) . we indeed recover the values in (lA.Gp . 

Finally, let us shortly discuss the non-Gaussian character of p Tz . comparing the 
expectation values Tr(p T 2 a m a n ), where m , n= 1,..., 4, with Tr(p T 2 aia 2 a 3 a 4 ), one observes 
that the Wick expansion, Eq. (T5]l . does not hold, unless c 2 = u\ v 2 . This remains true 
whatever basis we choose for the partial transpose. Thus, the partial transpose of a 
fermionic Gaussian state is usually not a Gaussian operator. 
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Here we consider the eigenvalue problem of the modified covariance matrices T±, that 
are associated to a pure-state covariance matrix F of the form ( I loll . The RDMs for 
subsystems A± and A 2 are determined via the reduced covariance matrices T 11 and 
T 22 . Since they split into two submatrices, one could equivalently consider the squared 
eigenvalue problem of the matrices 


= £ g m ,g„, = V M°fipm)rpn), ^(m) 

Ida p,q 


q (m) 

0 


m G A a , 
m £ A a , 


(B.l) 


with nonzero matrix elements only within m,n G A a , a 
M a have matrix elements 


m w = 5>p('W0, 

1 


m 

0 


1,2. The overlap matrices 


l G A a , 
l (f A a . 


(B.2) 


Note that both G 11 and G 22 have the same eigenvalues p\ < 1 with k = 

1 ,..., mindA!), \A 2 \), whereas p 2 k = 1 for the remaining eigenvalues of the larger matrix. 
We also introduce the block-diagonal matrix 


G 


G 11 0 \ 

0 G 22 J ’ 


(B.3) 


with all the nontrivial eigenvalues being doubly degenerate. 

The matrix elements of the covariance matrices r± in Eq. (1241) are determined 
through 

9mn = H, ( B - 4 ) 

with the vectors 


0> • (B.5) 


Thus the spectrum of Tj. follows from the eigenvalues of the squared matrix 

GL = ~ • ( B '6) 

P,Q 

Inserting (IB.5D and using the completeness property M^ q + M 2 q = 5 pq , the matrices G ± 
have the block form 


with 


G ± 


2 G 11 - 11 ±2 iF 

±2 iF T 2G 22 - 11 


(B.7) 


F m n = 


P,Q 


It is easy to check that F satisfies 
FF t = G n (]l- G 11 ) 


F F - G 22 (l - G 11 ), 


(B.8) 


(B.9) 
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and thus the following matrix identity holds 

(G ± -2G+ ll) 2 = -4G{1-G). (B.10) 

Rewriting in terms of the eigenvalues (n±) 2 and fi\ of G ± and G, respectively, one arrives 
at 

{VkY = 2 /4 - 1 ± 2 ¥h\j 1 - l4 - (^k Vlj ■ ( B - n ) 


Appendix C. Determinant formulas 

Let us consider the Gaussian operators 0± corresponding to the generalized covariance 
matrices r ± = tanh(W±/2). With a denoting the vector of Majorana operators, we 
introduce 


O (W ± ) = exp , 


Z(w±) = Trexp , 


(C.l) 


and thus 0± = 0(W±) /Z(W±) Our aim is to calculate various traces of the form 
Tr(0+O") with some integers m and n. Following the lines of Ref. [38], we first introduce 
the notation 


{W u W 2 } = Tr [O (Wi) O (W 2 )}. 


(C.2) 


We also note the simple fact that O m (W±) = O (mW±). Hence, the traces we consider 
can be written in the form 


TV(CT) = 


Z(mW + 


Tv (Of 1 Of) = 


{mW + , nW-} 


z m fw+y z m (w + )z n (w_)' 

They can be evaluated in terms of determinant formulas 


Z(W) = (±) t / det 2 cosh 


W 


{W U W 2 } = (±) v 7 det(l + e^'e^ 2 ), 


(C.3) 


(C.4) 


where the ± in parentheses symbolise the eventual sign ambiguity. The square root 
(and hence the sign ambiguity) indicates that the pairs of eigenvalues of the skew- 
symmetric matrices must be taken into account with halved degeneracy [35] . Note 
that, for Gaussian states commuting with the particle number operator (i.e., when the 
exponent can be written with a Hermitian matrix in terms of the fermion operators 
instead of Majoranas), similar trace formulas apply without square roots and sign 
ambiguity [58] . 

In Section [5] we need the traces of operators and O+O-, respectively. Applying 
Eqs. (IG.3D and (lG.4jh using hyperbolic identities for multiple arguments, one observes 
that the formulas can be expressed solely in terms of T-t- with the result 


Tr(O^) = ±4/det 


1 + 3 rl 


Tr (0 2 + 0_) = J det 


i + r 2 + 2T + r_ 


(C.5) 


The sign ambiguity can be fixed by comparing to exact calculations of the traces. We find 
that the negative sign in the first trace of Eq. (IG.5D is needed only in case T + contains a 
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quadruplet of purely imaginary eigenvalues. For the Ising chain, this can happen only in 
the symmetry-broken phase, h < 1. The numerics for small chains shows that, gradually 
decreasing the value of h, the appearance of this quadruplet exactly coincides with the 
vanishing of the first determinant. Interestingly, the second determinant in Eq. (1C.51) 
always remains positive and thus no sign ambiguity appears there. 
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